Stata：空间权重矩阵的构造

作者：潘星宇 (清华大学)      
连享会：(知乎 | 简书 | 码云 | CSDN)  

2019暑期“实证研究方法与经典论文”专题班

特别说明

文中包含的链接在微信中无法生效。请点击本文底部左下角的【阅读原文】，转入本文【知乎版】。

连享会空间计量专题推文：

Stata15 版本新添加了空间计量的官方命令， 本文将结合一些操作实例带领大家系统地了解 stata15 的官方空间计量命令—— sp 家族命令。

1. 命令概览

Stata15 提供的空间计量命令主要基于空间自回归模型展开，主要有以下几类：

A. 空间数据准备

• spshape2dta  将 shapefile 文件转换为 Stata 格式

• spset  将现有数据声明为空间数据

• spbalance  使面板数据具有很强的平衡性

• spcompress  压缩状态格式的 shapefile 文件

B. 空间权重矩阵的构建和空间滞后项的生成

• spmatrix  创建、操作和导入/导出权重矩阵

• spgenerate  生成空间滞后变量(wx)

C. 空间计量的回归命令

• spregress 截面 SAR 模型

• spivregress  工具变量 SAR 模型

• spxtregress  面板 SAR 模型

D. 空间计量的检验命令

• estat moran 基于moran指数的检验

• spregress postestimation 截面 SAR 模型的检验命令

• spivregress postestimation 工具变量 SAR 模型的检验命令

• spxtregress postestimation 面板 SAR 模型的检验命令

2. 准备工作：空间数据准备与空间权重矩阵的构建

让我们从一个简单的例子开始，研究二氧化碳排放以及与贫困的关系。有关二氧化碳排放的官方统计数据可在线获取 ; 我们仅为此保留了 2014 年的数据（并删除了给出区域总数的行）以生成此文件 LA-CO2-2014.dta。

我们还将使用多维贫困指数（IMD），这是一种结合了英国人口普查的信息的贫困数据。英国的 shapefile 文件可以从这里下载（未包含北爱尔兰）。Shapefile 通常有一个 .zip 压缩文件，其中包含几个文件组成部分。对于 Shapefile 文件的介绍可以参考空间权重矩阵的构建。

下面我们开始来进行数据的导入和声明工作。

1. /* 数据的导入和声明 */

2. *- 解压缩 Shapefile 文件

3. unzipfile "Local_Authority_Districts_December_2016_Full_Clipped_Boundaries_in_Great_Britain.zip" // 解压缩文件内容

4. *- 数据的导入

5. spshape2dta Local_Authority_Districts_December_2016_Full_Clipped_Boundaries_in_Great_Britain

6. *- 数据的声明

7. use Local_Authority_Districts_December_2016_Full_Clipped_Boundaries_in_Great_Britain.dta , clear

8. encode lad16cd, generate(lacode)

9. spset lacode, modify replace

在上述过程中，Stata 通过 spshape2dta 命令读取 shapefile 数据并构造了两个 .dta 文件。

其中一个更大文件名为 ……_shp.dta，其中包含了边界和位置的详细信息。
另一个文件名 为 …….dta，包含了 shapefile 数据本身带有的行政区划编码 "lad16cd" 。我们打开不含 _shp 后缀的文件，将 "lad16cd" 编码转换为数值型lacode的数值变量，然后使用 spset 命令将其声明为空间数据。这一步非常重要，因为声明为空间数据后，我们就可以进一步构建空间权重矩阵乃至进行空间回归的操作了。

然后我们将二氧化碳排放数据和多维贫困数据按照代码进行空间匹配

1. /* 数据匹配 */

2. merge 1:1 lad16cd using "LA-CO2-2014.dta" // 匹配二氧化碳排放数据

3. keep if _merge==3

4. drop _merge

5. save "CO2-merged.dta", replace //这一步我们就完成了给二氧化碳排放数据的编码

6. import excel "File_10_ID2015_Local_Authority_District_Summaries.xlsx", ///

7. sheet("IMD") firstrow clear

8. keep LocalAuthorityDistrictcode2 IMDAveragescore

9. rename LocalAuthorityDistrictcode2 lad16cd

10. merge 1:1 lad16cd using "CO2-merged.dta" //匹配多维贫困数据

11. keep if _merge==3

12. drop _merge

接下来我们来进行空间权重矩阵的构建。spmatrix 命令是一个 Stata 官方提供的比较好用的构建空间权重矩阵的命令。关于空间权重矩阵构建的其他知识可以参考空间权重矩阵的构建。

我们有两个权重矩阵构建的选择：第一个是假设与接壤的辖区可以相互影响。这有时称为adjacency matrix 或者 contiguity matrix （邻接矩阵）。第二种选择是使相关性与反距离成比例，它使用shapefile变量_CX和_CY来计算局部权限中点之间的距离。我们将尝试两者并进行比较。

1. /* 构建空间权重矩阵 */

2. spmatrix create contiguity W, rook

3. spmatrix create idistance W2

然后我们就可以分别得到一个边相邻的邻接矩阵（命名为 W ）和反距离矩阵（命名为 W2）。

3. 截面空间回归：回归与图示

准备工作完成后，我们就可以进行空间回归工作了。需要注意的是，stata15 官方命令只支持 SAR, 即
空间滞后模型 因此需要更多模型拓展的小伙伴还是需要其他外部命令的协助。

我们首先用最简单的 OLS 模型进行基准回归用以对比:

1. /* 基准回归 */

2. generate logdom = log(domestictotal)

3. generate logpop = log(population)

4. generate pc = logdom-logpop // 生成人均二氧化碳排放量的对数指标

5. regress pc IMDAveragescore // 基准回归

6. predict regres, residuals

多维贫困指数与人均二氧化碳排放量的系数为 -0.0082（95％置信区间为 -0.0066 到 -0.0097，p <0.001）：贫困指数越高的地区的排放量越低。

在这里我们可以用残差的可视化图来反映空间模型的适用性问题：

1. /* 残差的可视化 */

2. scatter regres domest, msymbol(Oh) name(regres, replace)

3. grmap regres, clnumber(9) fcolor(BuYlRd) // 在地图上绘制残差

4. name(regresmap, replace)

图中淡黄色表示残差为零（完全拟合），深蓝色表示较大的负残差（模型高估数据）和较暗红色表示大的正残差（模型低估了数据）。从图中可以看出残差在空间上是聚类的。有一些红色和蓝色的连片的大区域（即高值-低值聚集的区域）。伦敦周边的蓝色地区表明，尽管一些当局的贫困程度相当高，而另一些贫困程度较低，但人均二氧化碳排放总体上低于拟合值，而在英格兰北部边缘则相反，人均二氧化碳排放总体要高于拟合值。

我们再采用 spregress 命令进行空间滞后模型的回归。若附加 gs2sls 选项，则执行广义两阶段最小二乘法，语法格式如下：

1. spregress depvar [indepvars] [if] [in], gs2sls [gs2sls_options]

各个选项的含义如下：

• dvarlag(spmatname)  选项为必要选项，表示要调用的空间权重矩阵

• ivarlag(spmatname: varlist)  表示采用的工具变量

• errorlag(spmatname) 表示误差项也具有空间相关性（可以理解为SEM模型的选项）

我们亦可附加 ml 选项以便执行 MLE 估计，具体语法如下：

1. spregress depvar [indepvars] [if] [in], ml [ml_options]

具体估计结果如下：

1. *-空间 SAR 模型回归（相邻矩阵）

2. spregress pc IMDAveragescore, gs2sls dvarlag(W)

3. Spatial autoregressive model Number of obs = 326

4. GS2SLS estimates Wald chi2(2) = 154.82

5. Prob > chi2 = 0.0000

6. Pseudo R2 = 0.2928

7. -----------------------------------------------------------------------------

8. pc | Coef. Std. Err. z P>|z| [95% Conf. Interval]

9. ----------------+------------------------------------------------------------

10. IMDAveragescore | -.0062574 .0008112 -7.71 0.000 -.0078472 -.0046675

11. _cons | .5901449 .0273667 21.56 0.000 .5365072 .6437826

12. ----------------+------------------------------------------------------------

13. W

14. ----------------+------------------------------------------------------------

15. pc | .1716607 .0343417 5.00 0.000 .1043522 .2389692

16. -----------------------------------------------------------------------------

17. Wald test of spatial terms: chi2(1) = 24.99 Prob > chi2 = 0.0000

在相邻空间权重矩阵下，多维贫困指数与人均二氧化碳排放量的相关系数为 -0.0063 （ 95％ 置信区间 -0.0047至 -0.0078，z = -7.7，伪R2 = 0.29）;同样，贫困指数越高的地区的排放量越低。其空间自相关系数为 0.1716607，说明人均二氧化碳排放与相邻地区的排放量呈现出正相关关系。

1. *-空间 SAR 模型回归（反距离矩阵）

2. spregress pc IMDAveragescore, gs2sls dvarlag(W)

3. Spatial autoregressive model Number of obs = 326

4. GS2SLS estimates Wald chi2(2) = 189.68

5. Prob > chi2 = 0.0000

6. Pseudo R2 = 0.3868

7. -------------------------------------------------------------------------------

8. pc | Coef. Std. Err. z P>|z| [95% Conf. Interval]

9. ----------------+--------------------------------------------------------------

10. IMDAveragescore | -.0084929 .0007202 -11.79 0.000 -.0099046 -.0070813

11. _cons | .8576374 .0248013 34.58 0.000 .8090278 .906247

12. ----------------+--------------------------------------------------------------

13. W2

14. ----------------+--------------------------------------------------------------

15. pc | -.3431917 .0439964 -7.80 0.000 -.429423 -.2569604

16. -------------------------------------------------------------------------------

17. Wald test of spatial terms: chi2(1) = 60.85 Prob > chi2 = 0.0000

在反距离空间权重矩阵下，多维贫困指数与人均二氧化碳排放量的相关系数为 -0.0085 （ 95％ 置信区间 -0.0071 至 - 0.0099，z = -11.7，伪 -R2 = 0.37 ）。

1. *-空间SAR模型的残差可视化

2. // 相邻矩阵的残差可视化

3. spregress pc IMDAveragescore, gs2sls dvarlag(W)

4. predict spres, residuals

5. scatter spres domest, msymbol(Oh) name(spres, replace)

6. grmap spres, clnumber(9) fcolor(BuYlRd) name(spresmap, replace) title("SAR (adjacency) residuals")

7. twoway (scatter spres regres, ms(Oh) ///

8. ytitle("SAR (adjacency) residuals") ///

9. xtitle("Linear regression residuals")) ///

10. (function y=x, range(regres)), ///

11. legend(off) name(resvres, replace)

12. // 反距离矩阵的残差可视化

13. spregress pc IMDAveragescore, gs2sls dvarlag(W2)

14. predict spres2, residuals

15. scatter spres2 domest, msymbol(Oh) name(spres2, replace)

16. grmap spres2, clnumber(9) fcolor(BuYlRd) ///

17. name(spresmap2, replace) ///

18. title("SAR (inverse-distance) residuals")

19. twoway (scatter spres2 regres, ms(Oh) ///

20. ytitle("SAR (inverse-distance) residuals") ///

21. xtitle("Linear regression residuals")) ///

22. (function y=x, range(regres)), ///

23. legend(off) name(resvres2, replace)

相邻空间权重矩阵的残差分布

反距离空间权重矩阵的残差分布

哪种类型的相关矩阵是更合意的？伪 R2 表明反距离矩阵更适合数据，但答案实际上取决于我们对数据上下文的理解以及我们试图回答的问题。英国的地方行政区的规模大不相同，因此两个行政区可能会被另一个行政区分开，从而不连续，但他们在地理上非常接近，这一点特别是在城市地区非常明显。而面积非常大，人口却比较稀少的农村地区可能有相反的差异：地域上相连但中心之间的距离很远。从图中我们可以看到反距离矩阵更好地适用于那些邻近矩阵模型的红色和低估的大型农村地区，以及伦敦及周边地区。

4. 截面空间回归：检验命令

这一节主要介绍两个命令：estat impact以及 lrtest ，前者用于对空间效应的分解，后者为似然比检验，检验模型的选择。

效应分解：

• 平均直接影响: 衡量每个观测单位自变量对于因变量的平均总影响

• 平均间接影响: 衡量某个观测单位对所有其他观测单位的影响

1. *-效应分解

2. spregress pc IMDAveragescore, gs2sls dvarlag(W)

3. estat impact

4. ------------------------------------------------------------------------------

5. | Delta-Met

6. | dy/dx Std. Err. z P>|z| [95% Conf. Interval]

7. ----------------+-------------------------------------------------------------

8. direct

9. IMDAveragescore | -.0062816 .0008096 -7.76 0.000 -.0078685 -.0046947

10. ----------------+-------------------------------------------------------------

11. indirect

12. IMDAveragescore | -.0009863 .0002044 -4.83 0.000 -.0013868 -.0005857

13. ----------------+-------------------------------------------------------------

14. total

15. IMDAveragescore | -.0072679 .000856 -8.49 0.000 -.0089457 -.0055901

16. ------------------------------------------------------------------------------

可以发现，在相邻空间权重矩阵下。直接效应为 -0.0062816，即贫困指数每上升1个单位，直接人均二氧化碳排放下降0.63%；间接效应为 -0.0009863，表示本地人均二氧化碳排放量的变化引起的周边地区的二氧化碳排放量变化。

模型的选择：我们在模型中加入

1. . spregress pc IMDAveragescore, ml errorlag(W) dvarlag(W)

2. . estimates store normalreg

3. . spregress pc IMDAveragescore, ml dvarlag(W)

4. . lrtest normalreg

5. Likelihood-ratio test LR chi2(1) = 59.46

6. (Assumption: normalreg nested in .) Prob > chi2 = 0.0000

这里显示加入空间滞后误差项的模型在小于0.01%的水平下显著。因此我们可以认为残差项也存在着非常显著的空间依赖性。因此加入空间滞后误差项是更为合意的。

5. 结合 spgen 命令模型拓展

在更多的情境下，我们也要考虑自变量X也存在着空间依赖的可能性。但 stata 官方并未提供 sdm
的相关命令，这时我们可以采用 spgen 命令生成空间权重矩阵 W 与自变量X的空间滞后项，再采用 spregress 命令进行回归，从而可以得到考虑了自变量空间依赖性的模型结果。需要注意的是：这一方法缺乏对一致性和渐近分布的数理推导，因此在理论上不够严谨，因此使用的时候需要慎重选择。

1. spshape2dta Local_Authority_Districts_December_2016_Full_Clipped_Boundaries_in_Great_Britain, replace

2. spmatrix create contiguity W, rook \\ 读取 shapefile文件并生成空间权重矩阵

3. spgenerate Wtransporttotal = W*transporttotal

4. Variable Obs Mean Std. Dev. Min Max

5. pc 326 0.545492 0.129761 0.024488 0.965792

6. IMDAverage~e 326 19.46362 8.004457 5.009 41.997

7. WIMDAverag~e 326 14.35601 7.139981 0 34.3326

8. spregress pc IMDAveragescore WIMDAveragescore,gs2sls dvarlag(W)

9. Spatial autoregressive model Number of obs = 326

10. GS2SLS estimates Wald chi2(3) = 141.17

11. Prob > chi2 = 0.0000

12. Pseudo R2 = 0.3017

13. ------------------------------------------------------------------------------

14. pc | Coef. Std. Err. z P>|z| [95% Conf. Interval]

15. -----------------+------------------------------------------------------------

16. pc |

17. IMDAveragescore | -.0080134 .000819 -9.78 0.000 -.0096186 -.0064083

18. WIMDAveragescore | .0028504 .0014383 1.98 0.048 .0000314 .0056695

19. _cons | .639564 .0229944 27.81 0.000 .5944958 .6846322

20. -----------------+------------------------------------------------------------

21. W |

22. pc | .0488928 .0526389 0.93 0.353 -.0542776 .1520632

23. ------------------------------------------------------------------------------

24. Wald test of spatial terms: chi2(1) = 0.86 Prob > chi2 = 0.3530

在反距离空间权重矩阵下，多维贫困指数与人均二氧化碳排放量的相关系数为 -0.0085 （95％ 置信区间 -0.0071 至 - 0.0099，z = -11.7，伪 -R2 = 0.37 ）。而多维贫困指数也呈现出一定的空间依赖性。但这里我们的空间Wald检验没有通过，因此加入WX的模型在这里并不适用。

6. 面板模型的推广

掌握了上述操作后，空间面板模型的操作就非常简单了，这里只需要同时进行空间数据的声明 spset 以及面板数据的声明 xtset。具体的操作命令如下：

1. *-空间面板操作命令

2. copy http://www.stata-press.com/data/r15/homicide_1960_1990.dta .

3. copy http://www.stata-press.com/data/r15/homicide_1960_1990_shp.dta .

4. use homicide_1960_1990

5. xtset _ID year //声明面板数据

6. spset //声明空间数据

7. spmatrix create contiguity W if year == 1990 //以1990年数据为基准构建空间权重矩阵

8. spxtregress hrate ln_population ln_pdensity gini i.year, re dvarlag(W) //空间滞后模型回归

具体的参数解释也与空间截面的SAR模型类似，同样地，也可以通过 spgenerate 命令加入自变量的空间滞后项进行空间回归。

参考文献

1. Spatial Analysis in Stata 15

2. 空间权重矩阵的构建

3. Stata15 空间计量手册

• help spregress

• help spgenerate

• help spxtregress

• help lrtest

关于我们

• 【Stata 连享会(公众号：StataChina)】由中山大学连玉君老师团队创办，旨在定期与大家分享 Stata 应用的各种经验和技巧。

• 公众号推文同步发布于 CSDN-Stata连享会 、简书-Stata连享会 和 知乎-连玉君Stata专栏。可以在上述网站中搜索关键词Stata或Stata连享会后关注我们。

• 点击推文底部【阅读原文】可以查看推文中的链接并下载相关资料。

• Stata连享会 精彩推文1  || 精彩推文2

联系我们

• 欢迎赐稿： 欢迎将您的文章或笔记投稿至Stata连享会(公众号: StataChina)，我们会保留您的署名；录用稿件达五篇以上，即可免费获得 Stata 现场培训 (初级或高级选其一) 资格。

• 意见和资料： 欢迎您的宝贵意见，您也可以来信索取推文中提及的程序和数据。

• 招募英才： 欢迎加入我们的团队，一起学习 Stata。合作编辑或撰写稿件五篇以上，即可免费获得 Stata 现场培训 (初级或高级选其一) 资格。

• 联系邮件： [e-mail]

往期精彩推文

---

欢迎加入Stata连享会(公众号: StataChina)